home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / bump.exe / BUMP.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-02-10  |  16.2 KB  |  753 lines
  1. /*               Function definitions for
  2.                        BUMP 1.1
  3.         the Beginners Understandable Matrix Package
  4.                     for Borland C++
  5. */
  6. #include <iostream.h>    // for cout
  7. #include <ctype.h>        // for isdigit()
  8. #include <stdio.h>        // for fgets()
  9. #include <alloc.h>        // for coreleft()
  10. #include <stdlib.h>        // for atof()
  11. #include <conio.h>        // for cprintf()
  12. #include <string.h>        // for strlen()
  13. #include "bump.h"
  14. extern unsigned _stklen = 64000;
  15. int vcount = 0;
  16. ///////////////// The Vector member functions  /////////////////////////
  17.  
  18. Vector::Vector(int anr, int anc, char atemp, int afr, int afc)
  19. {
  20.     if(anr != 1 && anc != 1){
  21.         cout << "Either the number of rows or the number of columns must be 1.\n";
  22.         return;
  23.         }
  24.     which = vcount++;
  25.     nr = anr;
  26.     nc = anc;
  27.     temp = atemp;
  28.     fr = afr;
  29.     fc = afc;
  30.     lr = fr + nr - 1;
  31.     lc = fc + nc - 1;
  32.     nelem = nr + nc - 1;
  33.     if(nelem < 0 || nelem > 16000){
  34.         cout << nelem << "is an illegal number of elements for a vector.\n";
  35.         return;
  36.         }
  37.     if(nelem == 0) 
  38.         v = 0;
  39.     else if((v = new float[nelem]) == 0)
  40.         cout << "Out of memory.\n";
  41.     }
  42. Vector::Vector(Vector& a)
  43. {
  44.     int i;
  45.  
  46.     which = vcount++;
  47.     nr = a.nr;
  48.     nc = a.nc;
  49.     temp = a.temp;
  50.     fr = a.fr;
  51.     fc = a.fc;
  52.     lr = a.lr;
  53.     lc = a.lc;
  54.     nelem = a.nelem;
  55.     /* Before the "operator " functions return, they must "destroy" the Temp
  56.     vector they created.  But since many are returning Temp, C++ calls this
  57.     function to make a copy of Temp.  That call is implicit; the programmer
  58.     does not explicitly do it.  After that call, the destructor is called,
  59.     also implicitly, to destroy the original Temp.  This copy will then be
  60.     returned.  The following device detects that the vector being copied is a
  61.     temporary vector doomed to immediate destruction, so it just captures the
  62.     v array as its own and sets a.v = 0, so that the destructor won't hurt
  63.     it.  */
  64.     
  65.     if(a.temp == 'y'){
  66.         v = a.v;
  67.         a.v = 0;
  68.         }
  69.     else{
  70.         if((v = new float[nelem]) == 0){
  71.             cout << "Out of memory trying to create vector" << which <<".\n";
  72.             exit(1);
  73.             }
  74.         for(i = 0; i < nelem; i++)
  75.             v[i] = a.v[i];
  76.         }
  77.     }
  78. Vector::~Vector()
  79. {
  80.     if(v != 0 ) delete v; 
  81.     }
  82. // Free heap memory
  83. void Vector::freeh()
  84. {
  85.     if(v != 0 ) delete v; 
  86.     v = 0;
  87.     }
  88. int Vector::ReadA(char *filename)
  89. {
  90.     FILE *fp;
  91.     int i,j,k;
  92.     char s[241],field[30];
  93.  
  94.     if((fp = fopen(filename,"rt")) == 0){
  95.         cerr << "Cannot open " << filename << ".\n";
  96.         return(ERR);
  97.         }
  98.     fgets(s,240,fp);
  99.     j = 0;
  100.     for(i = 0; i < nelem; i++){
  101.         if(getfloat(v[i],fp,s,j) == ERR){
  102.             cerr << "Out of data in " << filename << ".\n";
  103.             break;
  104.             }
  105.         }
  106.     fclose(fp);
  107.     return(OK);
  108.     }
  109. void Vector::Display(char *title,int width, int decimalPlaces)
  110. {
  111.     int i, nperline;
  112.     char fmt[20];
  113.  
  114.     if(strlen(title) > 0) printf("\n%s\n", title);
  115.  
  116.     nperline = 79/width;
  117.     sprintf(fmt,"%s%d.%df", "%",width,decimalPlaces);
  118.  
  119.     for(i = 0;i<nelem;i++){
  120.         printf(fmt,v[i]);
  121.         if(i % nperline == nperline - 1) printf("\n");
  122.         }
  123.     }
  124.  
  125. Vector& Vector::operator = (Vector& a)
  126. {
  127.     int i;
  128.     int s = (nelem < a.nelem) ? nelem : a.nelem;
  129.  
  130.     for(i = 0; i < s; i++){
  131.         v[i] = a.v[i];
  132.         }
  133.     a.freet();
  134.     return(*this);
  135.     }
  136. Vector Vector::operator ~ () // Transposition
  137. {
  138.     int i;
  139.     Vector Temp(nc,nr,'y',fc,fr);
  140.     for(i = 0; i < nelem; i++)
  141.         Temp.v[i] = v[i];
  142.     return(Temp);
  143.     }
  144. float& Vector::operator [](int i)
  145. {
  146.  
  147.     if(v==0){
  148.         cerr << "Error: Reference to a deleted vector.\n";
  149.         exit(1);
  150.         }
  151.     if(nc == 1){
  152.         if (fr > i || lr < i){
  153.             cerr << "Illegal vector index: " << i << "\n";
  154.             exit(1);
  155.             }
  156.         return(v[i - fr]);
  157.         }
  158.     else{
  159.         if (fc > i || lc < i){
  160.             cerr << "illegal vector index: " << i << "\n";
  161.             exit(1);
  162.             }
  163.         return(v[i -fc]);
  164.         }
  165.     }
  166. Vector Vector::operator - (Vector &a)
  167. {
  168.     int i;
  169.     if(fr != a.fr || lr != a.lr){
  170.         cout << "vector dimensions do not match in + operator.\n";
  171.         }
  172.  
  173.     Vector Temp(nr,nc,'y');
  174.     for(i = 0; i < nelem; i++)
  175.         Temp.v[i] = v[i] - a.v[i];
  176.     freet();
  177.     return (Temp);
  178.     }
  179. //////////////// Friends of Vector (only) /////////////////////////
  180. /* operator + and operator - are written differently to illustrate
  181.    the two methods.
  182. */
  183. Vector operator + (Vector &a, Vector &b)
  184. {
  185.     int i;
  186.     if(b.fr != a.fr || b.lr != a.lr){
  187.         cout << "vector dimensions do not match in + operator.\n";
  188.         }
  189.  
  190.     Vector Temp(a.nr,a.nc,'y');
  191.     for(i = 0; i < a.nelem; i++)
  192.         Temp.v[i] = a.v[i] + b.v[i];
  193.     a.freet();
  194.     b.freet();
  195.     return (Temp);
  196.     }
  197. Vector operator * (float a, Vector &b)
  198. {
  199.     int i;
  200.     
  201.     Vector Temp(b.nr,b.nc,'y');
  202.     for(i = 0; i < b.nelem; i++)
  203.         Temp.v[i] = a * b.v[i];
  204.     b.freet();
  205.     return (Temp);
  206.     }
  207. Vector operator / (Vector &a, float b)
  208. {
  209.     int i;
  210.     float recip;
  211.     
  212.     Vector Temp(a.nr,a.nc,'y');
  213.     recip = 0;
  214.     if(b == 0)    printf("Division by zero in Vector operator /.\n");
  215.     else recip = 1./b;
  216.     for(i = 0; i < a.nelem; i++)
  217.         Temp.v[i] = a.v[i]*recip;
  218.  
  219.     a.freet();
  220.     return (Temp);
  221.     }
  222. //////////////////////  THE MATRIX MEMBER FUNCTIONS /////////////////////////
  223. Matrix::Matrix(int anr, int anc, char atemp, int afr, int afc)
  224. {
  225.     int i;
  226.  
  227.     which = vcount++;
  228.     nr = anr;
  229.     nc = anc;
  230.     temp = atemp;
  231.     fr = afr;
  232.     fc = afc;
  233.     lr = fr + nr - 1;
  234.     lc = fc + nc - 1;
  235.     if((m = new float* [nr]) == 0) goto gripe;
  236.     m -= fr;
  237.     for(i = fr; i <= lr; i++){
  238.         if((m[i] = new float [nc]) == 0) goto gripe;
  239.         m[i] -= fc;
  240.         }
  241.     return;
  242.     gripe:
  243.         cout << "Out of memory on array number " << which << ".\n";
  244.         exit(1);
  245.     }
  246. Matrix::Matrix(Matrix& a)
  247. {
  248.     int i,j;
  249.     which = vcount++;
  250.     nr = a.nr;
  251.     nc = a.nc;
  252.     temp = a.temp;
  253.     fr = a.fr;
  254.     fc = a.fc;
  255.     lr = a.lr;
  256.     lc = a.lc;
  257.     // See the note at this point in Vector::Vector(Vector& a), which
  258.     // applies here also. 
  259.     if(a.temp == 'y'){
  260.         m = a.m;
  261.         a.m = 0;
  262.         }
  263.     else{
  264.         if((m = new float* [nr]) == 0) goto gripe;
  265.         m -= fr;
  266.         for(i = fr; i <= lr; i++){
  267.             if((m[i] = new float [nc]) == 0) goto gripe;
  268.             m[i] -= fc;
  269.             }
  270.         for(i = fr; i <= lr; i++){
  271.             for(j = fc; j <= lc; j++)
  272.                 m[i][j] = a.m[i][j];
  273.             }
  274.         }
  275.     return;
  276.     gripe:
  277.         cout << "Out of memory on array number " << which << ".\n";
  278.         exit(1);
  279.     }
  280. Matrix::~Matrix()
  281. {
  282.     freeh();
  283.     }
  284. void Matrix::freeh()
  285. {
  286.     int i;
  287.  
  288.     if(m == 0) return;
  289.     for(i = fr; i <= lr; i++)
  290.         delete m[i];
  291.     delete m;
  292.     m = 0;
  293.     }
  294.  
  295. int Matrix::ReadA(char *filename)
  296. {
  297.     FILE *fp;
  298.     int i,j,k;
  299.     char s[241];
  300.  
  301.     if((fp = fopen(filename,"rt")) == 0){
  302.         cerr << "Cannot open " << filename << ".\n";
  303.         return(ERR);
  304.         }
  305.     fgets(s,240,fp);
  306.     k = 0;
  307.     for(i = fr; i <= lr; i++){
  308.         for(j = fc; j <= lc; j++){
  309.             if(getfloat(m[i][j],fp,s,k) == ERR)    break;
  310.             }
  311.         }
  312.     fclose(fp);
  313.     return(OK);
  314.     }
  315. void Matrix::Display(char *title,int width, int decimalPlaces)
  316. {
  317.     int i, j, nperline;
  318.     char fmt[20];
  319.  
  320.     if(strlen(title) > 0) printf("\n%s\n", title);
  321.  
  322.     nperline = 79/width;
  323.     sprintf(fmt,"%s%d.%df", "%",width,decimalPlaces);
  324.     for(i = fr; i <= lr; i++){
  325.         for(j = fc; j <= lc; j++){    
  326.             printf(fmt,m[i][j]);
  327.             if((j - fc % nperline == nperline -1) && j < lc)
  328.                 printf("\n");
  329.             }
  330.         printf("\n");
  331.         }
  332.     }
  333. Matrix& Matrix::operator = (Matrix& a)
  334. {
  335.     int i,j,nrmin,ncmin;
  336.     nrmin = nr < a.nr ? nr : a.nr;
  337.     ncmin = nc < a.nc ? nc : a.nc;
  338.  
  339.     if(nr != a.nr || nc != a.nc)
  340.         cout << "Equating matrices of unequal size.\n";
  341.     for(i = 0; i < nrmin; i++){
  342.         for(j = 0; j < ncmin; j++)
  343.             m[fr+i][fc+j] = a.m[a.fr+i][a.fc+j];
  344.         }
  345.     a.freet();
  346.     return(*this);
  347.     }
  348. Matrix Matrix::operator ~ () //Transposition
  349. {
  350.     int i,j;
  351.  
  352.     Matrix Temp(nc,nr,'y',fc,fr);
  353.     for(i = fr; i <= lr; i++){
  354.         for(j = fc; j <= lc; j++)
  355.             Temp.m[j][i] = m[i][j];
  356.         }
  357.     return(Temp);
  358.     }        
  359. float * Matrix::operator [](int i)
  360. {
  361.     if(m==0){
  362.         cerr << "Error: Reference to a deleted matrix.\n";
  363.         exit(1);
  364.         }
  365.  
  366.     if(i < fr || i > lr ){
  367.         cerr << "illegal matrix row index: " << i << "." << "\n";
  368.         cerr << "returning first row. \n";
  369.         return m[fr];
  370.         }
  371.     return(m[i]);
  372.     }
  373. Matrix Matrix::operator ! () // Inversion
  374. {
  375.     double determinant;
  376.  
  377.     if(temp == 'n'){
  378.         Matrix Temp(*this);
  379.         Temp.temp = 'y';
  380.         Temp.invert();
  381.         return(Temp);
  382.         }
  383.     else{ // If this matrix is temporary, turn it into its inverse.
  384.         invert();
  385.         return(*this);
  386.         }
  387.     }
  388. // Turn a matrix into its inverse. The value returned is the determinant.  
  389. double Matrix::invert(int startrow, int endrow)
  390. {
  391.     double pivot,determinant;
  392.     int i, j, k;
  393.  
  394.     if(startrow == 0) startrow = fr;
  395.     if(endrow   == 0) endrow   = lr;
  396.     determinant = 1.;
  397.     for(i = startrow; i <= endrow; i++){
  398.         if((pivot = m[i][i])== 0.){
  399.             printf("Zero diagonal; matrix probably singular. Returning partly inverted matrix.\n");
  400.             return(determinant);
  401.             }
  402.         determinant *= pivot;
  403.         m[i][i] = 1.;
  404.         for (j = fc; j <= lc; j++)
  405.             m[i][j] /= pivot;
  406.         for (k = fr; k <= lr; k++){
  407.             if(k == i) continue;
  408.             pivot = m[k][i];
  409.             m[k][i] = 0.;
  410.             for(j = fc; j <= lc; j++)
  411.                 m[k][j] -= pivot*m[i][j];
  412.             }
  413.         }
  414.     return(determinant);
  415.     }
  416. ///////////////FRIENDS OF MATRIX AND VECTOR  ////////////////////
  417. /*  M << b stores the vector b into the first column or row
  418.     of the Matrix M.  It seems impossible to use = for this
  419.     operation because the compiler insists that "operator ="
  420.     must be member, not merely a friend of the class. Thus it
  421.     is impossible to get at the protected elements of two different
  422.     classes in an "operator =" function.  Hence the use of <<.
  423.     */
  424. void operator << (Matrix& a, Vector& b)
  425. {
  426.     int i, nrmin, ncmin;
  427.     nrmin = min(a.nr,b.nr);
  428.     ncmin = min(a.nc,b.nc);
  429.     if(a.nr != b.nr || a.nc != b.nc)
  430.         cout << "<< applied to arrays of unequal size.\n";
  431.     if(b.nc == 1){  //The vector is a column
  432.         for(i=0;i<nrmin;i++)
  433.             a.m[a.fr+i][a.fc] = b.v[i];
  434.         }
  435.     else{
  436.         for(i= 0;i<ncmin;i++)
  437.             a.m[a.fr][a.fc+i] = b.v[i];
  438.         }
  439.     b.freet();
  440.     }
  441.  
  442. Matrix operator + (Matrix &a, Matrix &b)
  443. {
  444.     int i,j,fr,lr,fc,lc;
  445.     if(b.fr != a.fr || b.lr != a.lr){
  446.         cout << "matrix dimensions do not match in + operator.\n";
  447.         cout << "adding overlapping areas into a matrix with dimensions of the left one.\n";
  448.         }
  449.     fr = max(a.fr,b.fr);
  450.     lr = min(a.lr,b.lr);
  451.     fc = max(a.fc,b.fc);
  452.     lc = min(a.lc,b.lc);
  453.     
  454.     Matrix Temp(a.nr,a.nc,'y');
  455.     for(i = fr; i <= lr; i++){
  456.         for(j = fc; j <= lc; j++)
  457.             Temp.m[i][j] = a.m[i][j] + b.m[i][j];
  458.         }
  459.     a.freet();
  460.     b.freet();
  461.     return(Temp);
  462.     }
  463. Matrix operator - (Matrix &a, Matrix &b)
  464. {
  465.     int i,j,fr,lr,fc,lc;
  466.     if(b.fr != a.fr || b.lr != a.lr){
  467.         cout << "matrix dimensions do not match in + operator.\n";
  468.         cout << "adding overlapping areas into a matrix with dimensions of the left one.\n";
  469.         }
  470.     fr = max(a.fr,b.fr);
  471.     lr = min(a.lr,b.lr);
  472.     fc = max(a.fc,b.fc);
  473.     lc = min(a.lc,b.lc);
  474.     
  475.     Matrix Temp(a.nr,a.nc,'y');
  476.     for(i = fr; i <= lr; i++){
  477.         for(j = fc; j <= lc; j++)
  478.             Temp.m[i][j] = a.m[i][j] - b.m[i][j];
  479.         }
  480.     a.freet();
  481.     b.freet();
  482.     return(Temp);
  483.     }
  484. Matrix operator * (float a, Matrix& b)
  485. {
  486.     int i,j,fr,lr,fc,lc;
  487.     
  488.     Matrix Temp(b.nr,b.nc,'y'); 
  489.     for(i = b.fr; i <= b.lr; i++){
  490.         for(j = b.fc; j <= b.lc; j++)
  491.             Temp.m[i][j] = a*b.m[i][j];
  492.         }
  493.     b.freet();
  494.     return(Temp);
  495.     }
  496. Matrix operator / (Matrix& a, float b)
  497. {
  498.     int i,j,fr,lr,fc,lc;
  499.     float recip;
  500.  
  501.     recip = 0;
  502.     if(b != 0) recip = 1./b;
  503.     else {
  504.         printf("Dividing a matrix by 0. Result is 0.\n");
  505.         }
  506.     Matrix Temp(a.nr,a.nc,'y'); 
  507.     for(i = a.fr; i <= a.lr; i++){
  508.         for(j = a.fc; j <= a.lc; j++)
  509.             Temp.m[i][j] = recip*a.m[i][j];
  510.         }
  511.     a.freet();
  512.     return(Temp);
  513.     }
  514. Matrix operator * (Matrix &a, Matrix &b)
  515. {
  516.     int i,j,k,fr,lr,fc,lc;
  517.     double sum;
  518.  
  519.     Matrix Temp(a.nr,b.nc,'y',a.fr,b.fc); 
  520.     for(i = a.fr; i <= a.lr; i++){
  521.         for(j = b.fc; j <= b.lc; j++)
  522.             Temp.m[i][j] = 0.;
  523.         }
  524.     if(a.nc != b.nr){
  525.         cout << "matrix dimensions do not match in * operator.\n";
  526.         cout << "returning a zero matrix.\n";
  527.         }
  528.     else {
  529.         for(i = a.fr; i <= a.lr; i++){
  530.             for(j = b.fc; j <= b.lc; j++){
  531.                 sum = 0;
  532.                 for(k = 0; k < a.nc; k++)
  533.                     sum += a.m[i][a.fc+k]*b.m[b.fr+k][j];
  534.                 Temp.m[i][j] = sum;
  535.                 }
  536.             }
  537.         }
  538.     a.freet();
  539.     b.freet();
  540.     return(Temp);
  541.     }
  542. Vector operator * (Vector& a, Matrix& b)
  543. {
  544.     int i,j,fr,lr,fc,lc;
  545.     double sum;
  546.  
  547.     Vector Temp(a.nr,b.nc,'y',a.fr,b.fc); 
  548.     for(i = 0; i < Temp.nelem; i++){
  549.         Temp.v[i] = 0.;
  550.         }
  551.     if(a.nc != b.nr){
  552.         cout << "vector*matrix dimensions do not match.\n";
  553.         cout << "returning a zero vector.\n";
  554.         }
  555.     else {
  556.         for(i = 0; i < a.nelem; i++){
  557.             sum = 0;
  558.             for(j = 0; j < b.nr ; j++)
  559.                 sum += a.v[j]*b.m[b.fr+j][b.fc+i];
  560.             Temp.v[i] = sum;
  561.             }
  562.         }
  563.     a.freet();
  564.     b.freet();
  565.     return(Temp);
  566.     }
  567. Vector operator * (Matrix& a, Vector& b)
  568. {
  569.     int i,j,fr,lr,fc,lc;
  570.     double sum;
  571.  
  572.     Vector Temp(a.nr,b.nc,'y',a.fr,b.fc); 
  573.     for(i = 0; i < Temp.nelem; i++){
  574.         Temp.v[i] = 0.;
  575.         }
  576.     if(a.nc != b.nr){
  577.         cout << "matrix*vector dimensions do not match.\n";
  578.         cout << "returning a zero vector.\n";
  579.         }
  580.     else {
  581.         for(i = 0; i < Temp.nelem; i++){
  582.             sum = 0;
  583.             for(j = 0; j < a.nc; j++)
  584.                 sum += a.m[a.fr+i][a.fc+j]*b.v[j];
  585.             Temp.v[i] = sum;
  586.             }
  587.         }
  588.     a.freet();
  589.     b.freet();
  590.     return(Temp);
  591.     }
  592. Matrix operator * (Vector& a, Vector& b)
  593. {
  594.     int i,j,k,fr,lr,fc,lc;
  595.     double sum;
  596.  
  597.     Matrix Temp(a.nr,b.nc,'y',a.fr,b.fc); 
  598.     for(i = a.fr; i <= a.lr; i++){
  599.         for(j = b.fc; j <= b.lc; j++)
  600.             Temp.m[i][j] = 0.;
  601.         }
  602.     if(a.nc != b.nr){
  603.         cout << "vector*vector dimensions do not match.\n";
  604.         cout << "returning a zero matrix.\n";
  605.         }
  606.     else {
  607.         for(i = a.fr; i <= a.lr; i++){
  608.             for(j = b.fc; j <= b.lc; j++){
  609.                 sum = 0;
  610.                 for(k = 0; k < a.nc; k++)
  611.                     sum += a.v[k]*b.v[k];
  612.                 Temp.m[i][j] = sum;
  613.                 }
  614.             }
  615.         }
  616.     a.freet();
  617.     b.freet();
  618.     return(Temp);
  619.     }
  620. // A/B gives what is usually written as A'B, but without creating A'.
  621. Matrix operator / (Matrix& a, Matrix& b)
  622. {
  623.     int i,j,k,n;
  624.     double sum;
  625.  
  626.     Matrix Temp(a.nc,b.nc,'y',a.fc,b.fc);
  627.  
  628.     if(a.nr != b.nr){
  629.         cerr << "Matrices do not have same number of rows in / operator.\n";
  630.         cerr <<  "Using the smaller number.\n";
  631.         }
  632.     n = min(a.nr,b.nr);
  633.     for (i = a.fc; i <= a.lc; i++){
  634.         for(j = b.fc; j <= b.lc; j++){
  635.             sum = 0;
  636.             for(k = 0; k < n; k++)
  637.                 sum += a.m[a.fr+k][i]*b.m[b.fr+k][j];
  638.             Temp.m[i][j] = sum;
  639.             }
  640.         }
  641.     a.freet();
  642.     b.freet();
  643.     return(Temp);
  644.     }
  645. Vector operator / (Matrix& a, Vector& b)
  646. {
  647.     int i,j,k,n;
  648.     double sum;
  649.  
  650.     Vector Temp(a.nc,b.nc,'y',a.fc,b.fc);
  651.  
  652.     if(a.nr != b.nr){
  653.         cerr << "Matrices do not have same number of rows in / operator.\n";
  654.         cerr <<  "Using the smaller number.\n";
  655.         }
  656.     n = min(a.nr,b.nr);
  657.     for (i = a.fc; i <= a.lc; i++){
  658.         sum = 0;
  659.         for(k = 0; k < n; k++)
  660.             sum += a.m[a.fr+k][i]*b.v[k];
  661.         Temp.v[i-a.fc] = sum;
  662.         }
  663.     a.freet();
  664.     b.freet();
  665.     return(Temp);
  666.     }
  667. Vector operator / (Vector& a, Matrix& b)
  668. {
  669.     int i,j,k,n;
  670.     double sum;
  671.  
  672.     Vector Temp(a.nc,b.nc,'y',a.fc,b.fc);
  673.  
  674.     if(a.nr != b.nr){
  675.         cerr << "Vector and Matrix do not have same number of rows in / operator.\n";
  676.         cerr <<  "Using the smaller number.\n";
  677.         }
  678.     n = min(a.nr,b.nr);
  679.     for(j = b.fc; j <= b.lc; j++){
  680.         sum = 0;
  681.         for(k = 0; k < n; k++)
  682.             sum += a.v[k]*b.m[b.fr+k][j];
  683.         Temp.v[j-b.fc] = sum;
  684.         }
  685.     a.freet();
  686.     b.freet();
  687.     return(Temp);
  688.     }
  689. Matrix operator / (Vector& a, Vector& b)
  690. {
  691.     int i,j,k,n;
  692.     double sum;
  693.  
  694.     Matrix Temp(a.nc,b.nc,'y',a.fc,b.fc);
  695.  
  696.     if(a.nr != b.nr){
  697.         cerr << "Vectors do not have same number of rows in / operator.\n";
  698.         cerr <<  "Using the smaller number.\n";
  699.         }
  700.     n = min(a.nr,b.nr);
  701.     for (i = a.fc; i <= a.lc; i++){
  702.         for(j = b.fc; j <= b.lc; j++){
  703.             sum = 0;
  704.             for(k = 0; k < n; k++)
  705.                 sum += a.v[i-a.fc+k]*b.v[j-b.fc+k];
  706.             Temp.m[i][j] = sum;
  707.             }
  708.         }
  709.     a.freet();
  710.     b.freet();
  711.     return(Temp);
  712.     }
  713. ////////////  Utility Functions ///////////////////////////////////
  714. int getfloat(float& f, FILE *fp, char *s, int& j)
  715. {
  716.     int k;
  717.     char field[30];
  718.  
  719.     // skip over non-numerical characters
  720.     search: while(isnum(s[j]) == 0 && s[j] != '\n') j++;
  721.     if(s[j] == '\n'){
  722.         if(fgets(s,240,fp) == NULL){
  723.             return(ERR); // Hit end of file.
  724.             }
  725.         j = 0;
  726.         goto search;
  727.         }
  728.     k = 0;
  729.     while(isnum(s[j+k]) != 0 ){
  730.         field[k] = s[j+k];
  731.         k++;
  732.         }
  733.     field[k] = '\0';
  734.     j += k;
  735.     f = atof(field);
  736.     return(OK);
  737.     }
  738. // isnum returns 1 if c is  a digit, '.', '-', '+'.
  739. int isnum(char c)
  740. {
  741.     if (isdigit(c)) return 1;
  742.     if(c == '.' || c == '-' || c == '+' ) return 1;
  743.     return 0;
  744.     }
  745. int max(int a, int b)
  746. {
  747.     return(a < b ? b : a);
  748.     }
  749. int min(int a, int b)
  750. {
  751.     return(a < b ? a : b);
  752.     }
  753.